Multiple precision evaluation of the Airy Ai function with reduced cancellation

The series expansion at the origin of the Airy function Ai(x) is alternating and hence problematic to evaluate for x > 0 due to cancellation. Based on a method recently proposed by Gawronski, M\"uller, and Reinhard, we exhibit two functions F and G, both with nonnegative Taylor expansions at the origin, such that Ai(x) = G(x)/F(x). The sums are now well-conditioned, but the Taylor coefficients of G turn out to obey an ill-conditioned three-term recurrence. We use the classical Miller algorithm to overcome this issue. We bound all errors and our implementation allows an arbitrary and certified accuracy, that can be used, e.g., for providing correct rounding in arbitrary precision

The functions erf and erfc computed with arbitrary precision and explicit error bounds

The version available on the HAL server is slightly different from the published version because it contains full proofs.International audienceThe error function erf is a special function. It is widely used in statistical computations for instance, where it is also known as the standard normal cumulative probability. The complementary error function is defined as erfc(x)=erf(x)-1. In this paper, the computation of erf(x) and erfc(x) in arbitrary precision is detailed: our algorithms take as input a target precision t' and deliver approximate values of erf(x) or erfc(x) with a relative error bounded by 2^(-t'). We study three different algorithms for evaluating erf and erfc. These algorithms are completely detailed. In particular, the determination of the order of truncation, the analysis of roundoff errors and the way of choosing the working precision are presented. The scheme used for implementing erf and erfc and the proofs are expressed in a general setting, so they can directly be reused for the implementation of other functions. We implemented the three algorithms and studied experimentally what is the best algorithm to use in function of the point x and the target precision t'

A certified infinite norm for the implementation of elementary functions

The version available on HAL is slightly different from the published version because it contains full proofs.International audienceThe high-quality floating-point implementation of useful functions f : R -> R, such as exp, sin, erf requires bounding the error eps = (p-f)/f of an approximation p with regard to the function f. This involves bounding the infinite norm ||eps|| of the error function. Its value must not be underestimated when implementations must be safe. Previous approaches for computing infinite norm are shown to be either unsafe, not sufficiently tight or too tedious in manual work. We present a safe and self-validating algorithm for automatically upper- and lower-bounding infinite norms of error functions. The algorithm is based on enhanced interval arithmetic. It can overcome high cancellation and high condition number around points where the error function is defined only by continuous extension. The given algorithm is implemented in a software tool. It can generate a proof of correctness for each instance on which it is run

Évaluation efficace de fonctions numériques - Outils et exemples

With computers, it is possible to evaluate some numerical functions such as f = exp, sin, arccos, etc. The purpose of this thesis is to study how these functions can be implemented. Depending on the target architecture (software or hardware, small or high accuracy required), different problems must be addressed, but the final goal is always to eventually obtain an implementation as efficient as possible. We first study with an example the problems that arise in the case when the precision is arbitrary. When, on the contrary, the precision is known in advance, the function f is often replaced by an approximation polynomial p. Such a polynomial can then be very efficiently evaluated. In practice, the coefficients of p must be representable on a given finite number of bits. We propose several algorithms (some of them are heuristic and others are rigorous) for finding very good approximation polynomials satisfying this constraint. Our results also apply in the case when the approximant is a rational fraction. Once p has been found, one must prove that the error |p-f| is not greater than a given bound. The particular form of the function p-f (subtraction between two very close functions) makes this property hard to rigorously prove. We propose an algorithm for overcoming this difficulty. All these algorithms have been integrated into a software tool called Sollya, developed during the thesis. In the beginning, it was created for making the implementation of functions easier. Now, it may be interesting for anyone who needs to perform numerical computations in a safe environment.Les systèmes informatiques permettent d'évaluer des fonctions numériques telles que f = exp, sin, arccos, etc. Cette thèse s'intéresse au processus d'implémentation de ces fonctions. Suivant la cible visée (logiciel ou matériel, faible ou grande précision), les problèmes qui se posent sont différents, mais l'objectif est toujours d'obtenir l'implémentation la plus efficace possible. Nous étudions d'abord, à travers un exemple, les problèmes qui se posent dans le cas où la précision est arbitraire. Lorsque, à l'inverse, la précision est connue d'avance, la fonction f est souvent remplacée par un polynôme d'approximation p. Un tel polynôme peut ensuite être évalué très efficacement en machine. En pratique, les coefficients de p doivent être représentables sur un nombre fini donné de bits. Nous proposons un ensemble d'algorithmes (certains sont heuristiques, d'autres rigoureux) pour trouver de très bons polynômes d'approximation répondant à cette contrainte. Ces résultats s'étendent au cas où la fonction d'approximation est une fraction rationnelle. Une fois p trouvé, il faut prouver que l'erreur |p-f| n'excède pas un certain seuil. La nature particulière de la fonction p-f (soustraction de deux fonctions très proches) rend cette propriété difficile à prouver rigoureusement. Nous proposons un algorithme capable de contourner cette difficulté. Tous ces algorithmes ont été intégrés au logiciel Sollya, développé pendant la thèse. À l'origine conçu pour faciliter l'implémentation de fonctions, ce logiciel s'adresse à présent à toute personne souhaitant faire des calculs numériques dans un cadre complètement fiable

Evaluating a constant expression in multiple precision with a guaranteed error bound

The evaluation of special functions often involves the evaluation of numerical constants. When the precision of the evaluation is known in advance (e.g., when developing libms) these constants are simply precomputed once for a while. In contrast, when the precision is dynamically chosen by the user (e.g., in multiple precision libraries), the constants must be evaluated on the fly at the required precision and with a rigorous error bound. Often, such constants are easily obtained by means of formulas involving simple numbers and functions. In principle, it is not a difficult work to write multiple precision code for evaluating such formulas with a rigorous roundoff analysis: one only has to study how roundoff errors propagate through subexpressions. However, this work is painful and error-prone and it is difficult for a human being to be perfectly rigorous in this process. Moreover, the task quickly becomes impractical when the size of the formula grows. In this article, we present an algorithm that takes as input a constant formula and that automatically produces code for evaluating it in arbitrary precision with a rigorous error bound. It has been implemented in the Sollya free software tool and its behavior is illustrated on several examples.L'évaluation de fonction spéciales nécessite souvent d'évaluer certaines constantes. Lorsque la précision est connue à l'avance (par exemple, lorsqu'on développe une libm), ces constantes sont simplement précalculées une fois pour toutes. Mais lorsque la précision est fixée par l'utilisateur au moment de l'évaluation (comme c'est le cas pour les bibliothèques en précision arbitraire), les constantes doivent être évaluées à la volée à la précision demandée et en bornant rigoureusement les erreurs. Souvent, ce genre de constantes est donné par des formules faisant intervenir des fonctions simples et des entiers. En principe, écrire du code en précision arbitraire pour évaluer ce genre de formule avec une analyse rigoureuse des erreurs d'arrondi n'est pas une tâche difficile. Il suffit d'étudier comment les erreurs d'arrondi se propagent à travers les sous-expressions. Cependant, ce travail est ingrat, délicat, et il est difficile pour un être humain de rester parfaitement rigoureux. De plus, la tâche devient vite inabordable lorsque la formule grossit. Dans cet article, nous présentons un algorithme qui prend en entrée une formule constante et qui produit automatiquement du code pour l'évaluer en précision arbitraire et avec une borne d'erreur rigoureuse. Nous l'avons implémenté de façon expérimentale dans l'outil libre Sollya, et nous illustrons son comportement sur plusieurs exemples

Certified and fast computation of supremum norms of approximation errors

The version available on HAL corresponds to the version initially submitted to the conference and slightly differs from the published version since it does not account for remarks made by the referees.International audienceIn many numerical programs there is a need for a high-quality floating-point approximation of useful functions f, such as exp, sin, erf. In the actual implementation, the function is replaced by a polynomial p, leading to an approximation error (absolute or relative) epsilon = p-f or epsilon = p/f-1. The tight yet certain bounding of this error is an important step towards safe implementations. The main difficulty of this problem is due to the fact that this approximation error is very small and the difference p-f is highly cancellating. In consequence, previous approaches for computing the supremum norm in this degenerate case, have proven to be either unsafe, not sufficiently tight or too tedious in manual work. We present a safe and fast algorithm that computes a tight lower and upper bound for the supremum norms of approximation errors. The algorithm is based on a combination of several techniques, including enhanced interval arithmetic, automatic differentiation and isolation of the roots of a polynomial. We have implemented our algorithm and timings on several examples are given

Computation of the error functions erf & erfc in arbitrary precision with correct rounding

In this paper, the computation of erf(x) in arbitrary precision is detailed. A feature of our implementation is correct rounding: the returned result is the exact result (as if it were computed with infinite precision) rounded according to the specified rounding mode. The algorithm that computes the correctly rounded value of erf(x) for any argument x is detailed in this paper. In particular, the choice of the approximation formula, the determination of the truncation rank and of the computing precision are presented. When the current truncation rank and computing precision do not suffice to determine the correctly rounded value, they are increased and the computation is done again. The optimal strategy to adapt them is given in this paper. Finally, timings on some experiments are given, and the implementation of the complementary error function erfc is then sketched

Bilan carbone d'une équipe de recherche à l'Inria

Nous avons réalisé le bilan des émissions de gaz à effet de serre (GES) d’une équipe du centre de recherche Inria de Sophia Antipolis pour l’année 2019. L’équipe comprenait cette année-là 5 permanents, 5 doctorants, 1 post-doc, une assistante d’équipe à mi-temps et a accueilli pendant quelques mois deux stagiaires de master

Efficient and accurate computation of upper bounds of approximation errors

International audienceFor purposes of actual evaluation, mathematical functions f are commonly replaced by approximation polynomials p. Examples include floating-point implementations of elementary functions, quadrature or more theoretical proof work involving transcendental functions. Replacing f by p induces a relative error epsilon = p/f - 1. In order to ensure the validity of the use of p instead of f, the maximum error, i.e. the supremum norm of epsilon must be safely bounded above. Numerical algorithms for supremum norms are efficient but cannot offer the required safety. Previous validated approaches often require tedious manual intervention. If they are automated, they have several drawbacks, such as the lack of quality guarantees. In this article a novel, automated supremum norm algorithm with a priori quality is proposed. It focuses on the validation step and paves the way for formally certified supremum norms. Key elements are the use of intermediate approximation polynomials with bounded approximation error and a non-negativity test based on a sum-of-squares expression of polynomials. The new algorithm was implemented in the Sollya tool. The article includes experimental results on real-life examples

Dipole recovery from sparse measurements of its field on a cylindrical geometry

International audienceWe consider the inverse problem of recovering the position and moment of a magnetic dipolefrom sparse measurements of the field it generates, known on sections of three orthogonal cylindersenclosing it. This problem is motivated by recent measurements performed on Moon rocks, in viewof determining their magnetic properties. The key ingredient of the presented method is the useof rational approximation techniques, together with properties of the poles of the approximants, inorder to estimate the position of the dipole